# DR-white TIDE analysis script
#
# Purpose:
#   Compare Sanger trace files from an experimental DR-white sample with a
#   no-DSB parental DR-white control and estimate indel proportions.
#
# Required inputs:
#   - One control .ab1 or .scf file from parental DR-white without I-SceI
#   - One or more experimental .ab1 or .scf files
#   - The 20-bp target sequence downstream of the I-SceI recognition site
#
# User-editable fields near the end of this file:
#   Control    = control trace filename
#   ISceI      = 20-bp target sequence
#   samplename = one or more experimental trace filenames
#
# Required R packages:
#   Biostrings and sangerseqR (install from Bioconductor before running)
#
# Outputs:
#   - A graphical decomposition of each experimental trace
#   - A .csv file containing indel percentages and p values
#
# Interpretation for the canonical DR-white reporter:
#   - The 0-indel value represents no DSB or NHEJ without indels.
#   - The separate HR row represents the canonical 23-bp HR product.
#   - All other indel percentages are summed as NHEJ with processing.
#   - Total detectable repair = HR + NHEJ with processing.
#   This script is configured for the canonical DR-white reporter and has not
#   been validated for derivative reporters with different expected HR products.

import <- function(control_file, experimental_file, guide, seqstart = 100, seqend = 700, maxshift = 35, rg1 = NA, rg2 = NA) {
  ##Arguments: 
  ##control_file (char) & experimental_file (char): Sanger chromatogram files (.ab1 or .scf) 
  ##   control is typically DR-white with no I-SceI transgene
  ##   experimental has typically been treated with I-SceI transgene
  ## guide (char; internal argument name): 20 bp target sequence downstream of the 18 bp I-SceI recognition sequence in Sce.white
  ## seqstart (numeric): start of sequence read from where data will be included (because beginning of seq reads tends to be poor quality)
  ## seqend (numeric): last bp to be included in analysis (will be automatically adjusted if reads are shorter, see below)
  ## maxshift (numeric): range of basepair shifts (indels) to be analyzed, both positive and negative
  ## rg1, rg2 (numeric): [optional] the first (rg1) and last (rg2) base of the sequence region that is used for decomposition; will be automatically set if NA
  ##   Note: rg1&rg2 should be after the breaksite (if not, this will be corrected)
  
  require("Biostrings")
  require("sangerseqR")
  
  B<-c("A","C","G","T") #four bases, in the order that is always used by sangerseqR package
  
  patched.readsangerseq<-function(filename)
    #This is a slightly modified version of readsangerseq() in the sangerseqR package
    #It fixes a problem with reading some .ab1 files, which appear to have an aberrant last 
    #character in the sequence strings, which sangerseq() cannot cope with. It returns a sangerseq object.
  { require(sangerseqR)
    fc <- file(filename, open = "rb")
    rawdata <- readBin(fc, what = "raw", n = 1.2 * file.info(filename)$size)
    close(fc)
    filetype <- suppressWarnings(rawToChar(rawdata[1:4]))
    if (filetype == ".scf") {
      seq <- read.scf(filename)
    }
    else if (filetype == "ABIF") {
      seq <- read.abif(filename)
      l<-nchar(seq@data$PBAS.1)
      if(! substr(seq@data$PBAS.1,l,l) %in% LETTERS) { #if last character is not uppercase text
        seq@data$PBAS.1<-substr(seq@data$PBAS.1,1,l-1)
      }
      l<-nchar(seq@data$PBAS.2)
      if(! substr(seq@data$PBAS.2,l,l) %in% LETTERS) { #if last character is not uppercase text
        seq@data$PBAS.2<-substr(seq@data$PBAS.2,1,l-1)
      }
    }
    else stop("Invalid File.")
    return(sangerseq(seq))
  }
  
  ##load data (now automatically handles IBAF and SCF formats)
  control<-patched.readsangerseq(control_file)
  sample<-patched.readsangerseq(experimental_file)
  
  #extract primary sequences as called by sequencer:
  sequence_ctr <- primarySeq(control)
  sequence_mut <- primarySeq(sample)
  
  #adjust seqend to shortest sequence if necessary:
  seqend<-min(seqend,length(sequence_ctr), length(sequence_mut))
  
  
  #Alignments:
  
  #find position of I-SceI sequence (also if on opposite strand) and calculate breaksite: 
  Dguide<-DNAString(guide)
  if(!length(Dguide)==20){
    stop("target sequence should be 20 nucleotides")}
  guide.align.f <- pairwiseAlignment(pattern=Dguide, subject=sequence_ctr, type="local")
  guide.hit.f <- identical(as.character(Dguide), as.character(subject(guide.align.f))) #forward hit if full match found
  guide.align.r <- pairwiseAlignment(pattern=reverseComplement(Dguide), subject=sequence_ctr, type="local")
  guide.hit.r <- identical(as.character(reverseComplement(Dguide)), as.character(subject(guide.align.r))) #reverse hit if full match found  
  
  #there can only be one match with the top strands or with the bottom strand.
  if(guide.hit.f & !guide.hit.r) {breaksite<-start(subject(guide.align.f))+16}
  if(!guide.hit.f & guide.hit.r) {breaksite<-start(subject(guide.align.r))+3}
  if(guide.hit.f & guide.hit.r){stop("at least two target-sequence matches")}
  if(!guide.hit.f & !guide.hit.r){
    stop(paste("no target-sequence match 
               \n investigate whether the control sequence contains a Sanger sequencing mismatch. If so, change the target sequence to the corresponding IUPAC nucleotides in the control sequence.
               \n target sequence forward:", Dguide, 
               "\n target sequence reverse complement:", reverseComplement(Dguide),
               "\n control sequence:", sequence_ctr))}
  
  #align the sample to the control sequence and calculate offset
  if(seqstart>(breaksite-maxshift)){
              stop(paste("the breaksite (",breaksite,") is too close to the start of the sequence read -> If possible set start of sequence read lower.
              The sequence start of sequence read is maximal n bp breaksite - n bp of the chosen indel size range."))}
  
  #get the sequence interval on which the alignment should be based:
  if((nchar(sequence_ctr) < ((breaksite-maxshift)-seqstart)) | (nchar(sequence_mut) < ((breaksite-maxshift)-seqstart)))
  {
    stop(paste("(one of) the sequence run(s) is too short for proper alignment"))}
  
  seq_ctr <- substr(sequence_ctr, seqstart, breaksite-maxshift)
  seq_mut <- substr(sequence_mut, seqstart, breaksite-maxshift)
  align_seq <- pairwiseAlignment(pattern = seq_mut, subject = seq_ctr, type = "local")
  if(align_seq@score<20){
    stop("there is no good alignment found between the control amd test sample -> 
        The alignment window is too small or of bad quality or the control and test sample do not match")}
  offset_mut <- align_seq@pattern@range@start-align_seq@subject@range@start
    
  #extract control data:
  if (control@primarySeqID == "From scf file") {
    if(peakPosMatrix(control)[1,1]==0){peakPosMatrix(control)[1,1] <- 1}
    peak_ctr_loc <- peakPosMatrix(control)[,1]
  } else if (control@primarySeqID == "From ab1 file") { 
  if(peakPosMatrix(control)[1,1]==1){peakPosMatrix(control)[1,1] <- 2}
  peak_ctr_loc <- peakPosMatrix(control)[,1]-1 #for some reason sangerseq() added 1, so we substract it again
  }
  
  peak_ctr_height <- traceMatrix(control)[peak_ctr_loc,] #matrix with a column for each base 
  peak_ctr_height <- peak_ctr_height[1:seqend,]
  peak_ctr_height[is.na(peak_ctr_height)]<-0 #set NAs to 0
  colnames(peak_ctr_height)<-B
  
  #extract experimental data:
  if (sample@primarySeqID == "From scf file") {
    if(peakPosMatrix(sample)[1,1]==0){peakPosMatrix(sample)[1,1] <- 1}
    peak_mut_loc <- peakPosMatrix(sample)[,1]
  } else if (sample@primarySeqID == "From ab1 file") { 
    if(peakPosMatrix(sample)[1,1]==1){peakPosMatrix(sample)[1,1] <- 2}
  peak_mut_loc <- peakPosMatrix(sample)[,1]-1 #sangerseq() added 1, we substract it again
  }
  
  peak_mut_height <- traceMatrix(sample)[peak_mut_loc,] #matrix with a column for each base 
  peak_mut_height <- peak_mut_height[1:seqend,]
  peak_mut_height[is.na(peak_mut_height)]<-0 #set NAs to 0
  colnames(peak_mut_height)<-B
   
  #set rg1 and rg2, if not provided by user:
  rg1<-ifelse(is.na(rg1), breaksite+maxshift+5, rg1)
  rg2<-ifelse(is.na(rg2), seqend-maxshift-5, rg2)
  
  #check if rg1 and rg2 are within meaningful range: 
  if(rg1< breaksite+maxshift+5) {
    rg1<- breaksite+maxshift+5
    warning(paste("left boundary of decomposition window was adjusted", rg1, 
                  "It must be at least 5bp plus the maximum indel size downstream of the expected break site"))
  }
  
  if(rg2 > seqend-maxshift-5) {
    rg2<- seqend-maxshift-5
    warning(paste("right boundary of decomposition window was adjusted to",rg2,
                  "It cannot be more than the length of the shortest sequence read minus the maximum indel size minus 5bp."))
  }
  
  if(rg2 > seqend-offset_mut) {
    rg2<- seqend-offset_mut-5
    warning(paste("right boundary of decomposition window was adjusted to",rg2,
                  "It cannot be more than the length of the shortest sequence read minus the maximum indel size minus 5bp."))
  }
  
  if(rg2<rg1+maxshift*2) {
    stop(paste("boundaries of decomposed region are not acceptable -> 
        Set boundaries further apart or use smaller indel size if possible. 
        Maximum decomposition window spans from 5bp + n bp indel size range downstream of the break to 5bp + n bp indel size from the end of the shortest sequence read"))}
  
  #control for wrongly not/extra annotated peaks
  ctr_loc1 <- ctr_loc2 <- NA
  ctr_loc1<- peak_ctr_loc[1:seqend]
  ctr_loc2<- peak_ctr_loc[2:seqend]
  
  mut_loc1 <- mut_loc2 <- NA
  mut_loc1<- peak_mut_loc[1:seqend]
  mut_loc2<- peak_mut_loc[2:seqend]
  
  #average distance between peaks
  ctr_distance <- cbind(ctr_loc1[(seqstart+1):(seqend)]-ctr_loc1[(seqstart):(seqend-1)])
  ctr_av_dis <- colMeans(ctr_distance)
  
  mut_distance <- cbind(mut_loc1[(seqstart+1):(seqend)]-mut_loc1[(seqstart):(seqend-1)])
  mut_av_dis <- colMeans(mut_distance)
  
  #adjusted average distance for the smallest average distance
  if(ctr_av_dis<mut_av_dis) {
    Means_dis_s <- ctr_av_dis
  } else {
    Means_dis_s <- mut_av_dis
  }
  
  #adjusted average distance for the biggest average distance
  if(ctr_av_dis>mut_av_dis) {
    Means_dis_b <- ctr_av_dis
  } else {
    Means_dis_b <- mut_av_dis
  }
  
  #check for abnormalities in the entire sample, is the spacing between the nucleotides contant or not.
  ctr_outlier_s <- mut_outlier_s <- 0
  ctr_outlier_b <- mut_outlier_b <- 0
  
  ctr_outlier_s <- which(ctr_loc2[seqstart:rg2]>ctr_loc1[seqstart:rg2]+(3*Means_dis_s/2))+seqstart
  mut_outlier_s <- which(mut_loc2[seqstart:rg2]>mut_loc1[seqstart:rg2]+(3*Means_dis_s/2))+seqstart
  ctr_outlier_b <- which(ctr_loc2[seqstart:rg2]<ctr_loc1[seqstart:rg2]+(Means_dis_b/2))+seqstart
  mut_outlier_b <- which(mut_loc2[seqstart:rg2]<mut_loc1[seqstart:rg2]+(Means_dis_b/2))+seqstart
  
  if(length(ctr_outlier_s)>0 | length(mut_outlier_s)>0 | length(ctr_outlier_b)>0 | length(mut_outlier_b)>0){
    warning(paste("the spacing between the nucleotides in (one of) sanger sequence file(s) is not contant. This might indicate for wrongly not or extra annotated nucleotides. This can influence the TIDE estimation, check the chromotogram for abnormalities"))}
  
  
  return(list(
    ctr=peak_ctr_height, 
    mut=peak_mut_height, 
    seqstart=seqstart,
    seqend=seqend,
    maxshift=maxshift,
    rg1=rg1,
    rg2=rg2,
    breaksite=breaksite,
    offset_mut=offset_mut,
    B=B,
    experimental_file=experimental_file))
}

quality <- function(import) {
  ## All the arguments are generated in the function 'TIDE_import'. 
  ## import$ctr = peakheigths of the control sample (e.g. no I-SceI transgene)
  ## import$mut = peakheigths of the sample that have had a DSB/repair
  ## import$breaksite = site the I-SceI nuclease is supposed to break according to literature 
  ## import$seqstart = start of sequence read from where data will be included (because beginning of seq reads tends to be poor quality)
  ## import$seqend = last bp to be included in analysis (will be automatically adjusted if reads are shorter, see below)
  ## import$maxshift = which basepair shifts (indels) you want to know the percentage of.
  ##   Note: import$maxshift is the number to one direction, in the calculation it determines the shift to both direction (deletion & insertion) 
  ## import$rg1/import$rg2 = the sequence trace that is used for decomposition
  ##   Note: rg1&rg2 should be always after the breaksite
  # import$offset_mut = the offset that seuquence trace of the sample has with repect to the control sequence trace.
  
  ## The function will return a plot of the percentages of aberrant sequence trance per location.
  ##   plot will indicate expected breaksite location
  ##   plot will indicate the sequence window that is used for decomposition   
  ## The function will return the difference percentages of aberrant sequences compared to the control
  
  #Calculate percentage of each bp per peak and correct for the offset
  procent_ctr <- import$ctr  
  procent_mut <- import$mut;
  if (import$offset_mut>0){
    procent_mut <- rbind(import$mut[(1+import$offset_mut):nrow(import$mut),],matrix(NA,import$offset_mut,4))
    procent_mut <- (procent_mut/(rowSums(procent_mut)))*100
    procent_ctr <- (import$ctr/(rowSums(import$ctr)))*100
  } else if (import$offset_mut<0){
    procent_mut <- rbind(matrix(NA,-import$offset_mut,4),import$mut[1:(nrow(import$mut)+import$offset_mut),])
    procent_mut <- (procent_mut/(rowSums(procent_mut)))*100
    procent_ctr <- (import$ctr/(rowSums(import$ctr)))*100
  } else if (import$offset_mut==0){
    procent_ctr <- (import$ctr/(rowSums(import$ctr)))*100
    procent_mut <- (import$mut/(rowSums(import$mut)))*100
  }
  
  ## calculate total percentage mutations
  percentage_mutation_ctr <- rowSums(procent_ctr * t(apply(procent_ctr,1,function(x){!(x==max(x))})))
  percentage_mutation_sample <- rowSums(procent_mut * t(apply(procent_ctr,1,function(x){!(x==max(x))})))
  
  #plot aberrant sequence signal
  plot(percentage_mutation_sample, 
         type="h", col="green3", 
         xlim=c(import$seqstart, import$seqend), 
         ylim=c(0,100),
         xlab="basepair",
         ylab="% of aberrant sequences")
    lines(percentage_mutation_ctr, type="h", col="black")
    
    legend("topleft",legend=c("control sample", "test sample"), bty="n", pch=15, col=c(1,3))
    
    #show decomposition window
    rect(import$rg1, 110 , import$rg2, 110, density = 1, xpd=TRUE, col="grey", lwd=6) 
    text(import$rg1+((import$rg2-import$rg1)/2), (110+6), xpd=TRUE, col="grey", as.character("region for decomposition"))    
    
    #indicate theoretical breaksite
    if (import$breaksite>0){
      abline(v=import$breaksite, lty=5,lwd=3,col='blue')
      legend("topright",legend=paste('expected cut at ',import$breaksite,'bp',sep=''),text.col='blue', bty="n")}
    else{
      legend("topright",legend='no cut',text.col='blue', bty="n")
    }
   
  #calculate average mutation percentage 
  meanper_ctr_prebreak <- mean(percentage_mutation_ctr[import$seqstart:(import$breaksite-20)])
  meanper_mut_prebreak <- mean(percentage_mutation_sample[import$seqstart:(import$breaksite-20)])
  meanper_ctr_postbreak <- mean(percentage_mutation_ctr[import$breaksite:(import$seqend-20)])
  meanper_mut_postbreak <- mean(percentage_mutation_sample[import$breaksite:(import$seqend-20)]) 
  
  #print the percentages of each shift
  percentage_mutation <- data.frame(percentage = round(c(meanper_ctr_prebreak, meanper_mut_prebreak, meanper_ctr_postbreak, meanper_mut_postbreak),1))
  rownames(percentage_mutation) <- c("mean % pre-break control sample", "mean % pre-break test sample", "mean % post-break control sample", "mean % post-break test sample");

  print("");
  print(import$experimental_file);
  print(percentage_mutation)
  
}

decomposition <- function(import,  p.threshold = 0.001) {     
  ## All the arguments are generated in the function 'TIDE_import', except p.threshold
  ## import$ctr = peakheigths of the control sample (e.g. no I-SceI transgene)
  ## import$mut = peakheigths of the sample that have had a DSB/repair
  ## import$maxshift = size range of indels to be considered in the decomposition.
  ##   Note: calculation is always done in both directions, i.e. for both deletions and insertions of sizes 0:maxshift 
  ## import$rg1 = first base in the sequence sequence traces that is used for decomposition
  ## import$rg2 = last base in the sequence sequence traces that is used for decomposition
  ##   Note: rg1&rg2 should be always after the breaksite
  ## import$offset_mut = the offset that seuquence trace of the sample has with repect to the control sequence trace.
  ## p.threshold = p-value signicance threshold
  
  ## The function will generate a barplot with the prediction of the most prominent indels in the population of cells  
  ## The function will return the percentages of each indel in the sample with associated p-value
  
  require("colorspace")
  require("nnls")
  
  shiftrange<-c(-import$maxshift: import$maxshift) 
  
  #decomposite import$mut sequence data into indel combinations, 
  #separately for each base in c("A","C","G","T"). Stack up the data for the four bases in one aggragation matrix
  I_matrix <- c()
  I_vec <- c()
  
  for(b in import$B) #loop through four bases
  {#simulate sequencing peak data for all hypothesized indels from control peaks:
    sim <- matrix(NA, nrow=import$rg2-import$rg1+1, ncol=import$maxshift*2+1)
    colnames(sim)<-shiftrange
    for(i in shiftrange) {sim[,as.character(i)] <- import$ctr[(import$rg1:import$rg2)-i,b]}
    I_matrix <- rbind(I_matrix, sim)
    I_vec <- c(I_vec,import$mut[(import$rg1:import$rg2)+import$offset_mut,b])
  }
  
  #non-negative linear fit:
  NNFIT <- nnls(I_matrix,I_vec)
  
  ## pvalue calculation (source: https://www.princeton.edu/~slynch/soc504/mult_reg2.pdf)
  #standard error:
  se <- sqrt(diag((sum((NNFIT$fitted-I_vec)^2)/(nrow(I_matrix)-(import$maxshift*2+1)))*solve(t(I_matrix)%*%I_matrix)));
  
  #p-value:
  pv <- 2*pnorm(-abs(NNFIT$x/se))
  
  #R^2
  Rsq <- cor(NNFIT$fit,I_vec)^2
  
  #components in percentages:
  comper<-(Rsq*100*(NNFIT$x/sum(NNFIT$x)))

  #HR
  if (import$maxshift<23){print("Warning, this will probably crash because maxshift is less than my HR product size of 23")}
  exp_HR <- -23 #your HR result will lead to -23 deletion, we assume that all -23 deletions = HR event
  HRT <- which(shiftrange==exp_HR) 
  comper[import$maxshift*2+2]=comper[HRT]
  comper[HRT] <- 0
  pv[import$maxshift*2+2]=pv[HRT]
  pv[HRT] <- 1
   
  COL <- ifelse(pv<p.threshold,"red","black")
  COL[import$maxshift+1] <- ifelse(pv[import$maxshift+1]<p.threshold,"#FF000080","black")

  shiftnames <- c(shiftrange, import$maxshift+1)
  
  #plot decomposition graph              
  bp <- barplot(comper, 
                col=COL, 
                border = COL, 
                names.arg=shiftnames, 
                ylim=c(0, max(comper+10)), 
                xlab="<--deletion     insertion-->", 
                ylab="% of sequences",  
                xaxt='n')
    
    #make x-axis 
    a <- min(ceiling(shiftrange/5)*5)
    p <- pretty (c(a:-a), n=(round((length(shiftrange)-1)/5,0)-1))
    axis(1,at=bp[p+max(p)+1+a+max(shiftrange)],labels=p)
  
    axis(1,at=bp[max(shiftrange)*2+2.5],labels="HR")
    axis(1,at=bp[HRT],labels="NA")
    
    #above each group of bars: show percentage (mean across four bases)
  if(length(bp[pv<p.threshold]) > 0)
  {  
  text(bp[pv<p.threshold], (comper+5)[pv<p.threshold], as.character(((round(comper,1))[pv<p.threshold])))
  }  
    #display Rsq values as an indication of the accuracy:
    legend("topright",legend=as.expression(c(bquote(p < .(p.threshold)), bquote(p >= .(p.threshold)))), title= as.expression(bquote(R^2 == .(round(Rsq,2)))), pch=15, col=c("red",'black'), bty="n")
  
  eff <- round((Rsq*100) - comper[import$maxshift+1],1)
  cat("overall efficiency =", eff, "%\n")
  
  decomp.summary <- data.frame(percentage = round(comper,1), pvalue = signif(pv,2))
  rownames(decomp.summary) <- shiftnames
  rownames(decomp.summary)[which(shiftrange>0)]=paste('+',rownames(decomp.summary)[which(shiftrange>0)],sep='')
  
  rownames(decomp.summary)[length(shiftnames)]<- paste("HR")
  
  write.csv(file=paste(import$experimental_file, "decomp_summary.csv", sep ="_"),decomp.summary)
  
  invisible(list(
    bp=bp,
    pv=pv, 
    p.threshold=p.threshold, 
    NNFIT=NNFIT
    ))  
}

Control <- "control.ab1"
ISceI <- "CAAGATCCTTCTGATGGCCG"

samplename <- c("sample1.ab1", "sample2.ab1")
for (i in samplename){
	im <- import(control_file=Control, experimental_file=i, guide= ISceI, maxshift = 35,
	seqstart=50, rg1=300, rg2=400);
	q <- quality(im);
	d <- decomposition(im)
}